GRADIENT RECOVERY IN ADAPTIVE FINITE ELEMENT 
METHODS FOR PARABOLIC PROBLEMS 



OMAR LAKKIS AND TRISTAN PRYER 

Abstract. We derive energy-norm a posteriori error bounds, using gradient 
recovery (ZZ) estimators to control the spatial error, for fully discrete schemes 
for the linear heat equation. This appears to be the first completely rigorous 
derivation of ZZ estimators for fully discrete schemes for evolution problems, 
without any restrictive assumption on the timestep size. An essential tool for 
the analysis is the elliptic reconstruction technique. 

Our theoretical results are backed with extensive numerical experimenta- 
tion aimed at (a) testing the practical sharpness and asymptotic behaviour of 
the error estimator against the error, and (b) deriving an adaptive method 
based on our estimators. 

An extra novelty provided is an implementation of a coarsening error "preindi- 
cator", with a complete implementation guide in ALBERTA in the appendix. 



1. Introduction 

Gradient recovery a posteriori error estimators have been widely used since their 
dissemination in the engineering and scientific computation community by [ZZ87 , 
for which we will often refer to them shortly as ZZ estimators. Since their introduc- 
tion they have constituted the most serious rival to residual estimators introduced 
earlier on in |BR78| ■ The key to ZZ estimator's success is their implementation's 
simplicity, mild dependence upon the problem's data, and striking superconver- 
gence and asymptotic exactness properties. On the other hand, residual esti- 
mators, which are the main competitor to ZZ estimators, are a bit more involved 
in implementation and cost more to compute, but they are easier to handle from 
the mathematical analysis view-point in deriving rigorous upper and lower bounds. 
This situation has led to most of the theoretical results for evolution equations 
being obtained in the last two decades via residual estimators; we refer to [LM06 
for a review. Meanwhile rigorous mathematical work on recovery estimators has 
progressed, especially in the last decade, but mostly for stationary elliptic equa- 
tions, e.g., |AO00l lBX03al IBX03bl ICB02l IFV061 ILZ991 IPTcM IXZ04] . In c ontrast , 
very little progress was made on evolution problems, where an exception is [LW06 , 
where the main analytic difficulty comes from the singularly perturbed nature of 
the elliptic problems arising from time-stepping procedures. 

The aim of our work is to bridge the gap between the practical use of ZZ es- 
timators in adaptivity for evolution equations, studied by ZW98, Pic03J, and the 
rather mature error control theory via recovery for stationary equations. We fo- 
cus on the model problem provided by the linear heat equation. [LW06 , who are 
to our knowledge the only researchers to have explored this issue in depth, while 
obtaining satisfactory error bounds for spatially discrete schemes, must assume un- 
realistically small time-steps for the fully discrete case. In this paper we push one 
step forward by thoroughly analysing the fully discrete backward Euler schemes. 
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Namely, we provide reliable error bounds. The efficiency and asymptotic exactness 
of the bounds is dealt with computationally. 

Our main analytical tool to tackle the fully discrete scheme's difficulties is the el- 
liptic reconstruction in the fully discrete context, studied in [LM06 ], which provides 
a way to take advantage of elliptic a posteriori error estimates based on gradient 
recovery following the exposition of [AOOOj . 

The elliptic reconstruction technique, introduced under this name by |MN03j . 
involves the error's splitting into two parts, a parabolic error and an elliptic error, 
through the use of the elliptic reconstruction of the discrete solution, defined in 
(3.1). This allows to utilise existing elliptic a posteriori estimators for the elliptic 
part and standard parabolic energy estimates to control the second part. Despite 
this technique being initially introduced to derive sharp bounds for lower order 
spatial error norms, such as L 2 (i?) |MN03| ILM061 ILM07j and ^(0) |DLM09j . 
we apply it here as an analysis tool in an energy-norm framework, where a direct 
approach may lead to a highly complicated analysis for the fully discrete scheme. 
In fact, the single most interesting feature of the elliptic reconstruction, is that the 
parabolic error's energy norm term is of higher order (with respect to the spatial 
mesh-size parameter) than the elliptic error, as seen in [LM06J. In this paper we 
show, rigorously, that the full energy error can be accounted for only by the elliptic 
error, as long as data and time-step are resolved well enough (cf. Lemma 3.3). 
This crucial observation is also used to obtain residual a posteriori estimates for 
nonconforming methods in GL08 . Note that, it is part of the adaptive methods 
practitioner's folklore to employ heuristic versions of this argument. By way of ex- 
ample, we quote "the [full parabolic] discretisation in energy norm can be bounded 
by the [elliptic error] estimator" from [ZW98 . 

Although we treat the case of the Laplace operator, for simplicity, in this pa- 
per, our results can be extended to cover more general elliptic operators, even 
time-dependent ones, by using appropriate elliptic gradient recovery techniques, 
described by |FV06j . and a more careful time-step analysis, as in [GL08 . 

The paper is organised as follows. In <|2]we introduce the model problem, and its 
discretisations via conforming finite elements in space and backward Euler in time 
and we review the known results, about recovery estimators for elliptic problems, 
that will be used in the sequel. In f|3]we describe the elliptic reconstruction tech- 
nique and illustrate its use for the spatially semidiscrete problem. This paves the 
way to tackle the fully discrete problem in Sj4j where our main results are stated. 
In <|5 using numerical tests, we study the practical behaviour of the estimators and 
in S|6 we explore the adaptive schemes based on our estimators. 

As we have used the finite element toolbox ALBERTA, written and documented 
by |SS05j . for the tests, we have taken the opportunity to implement a coarsening 
preindicator, previously unavailable and (for space's sake) fully described in the 
Appendix [X] This estimator predicts the "information loss" error that will occur 
under coarsening of the mesh at each timestep of the adaptive method and is crucial 
in an adaptive code to control information loss during coarsening. 



2. Set up 

2.1. The model problem. Let fi C R d be a bounded polyhedral domain and 
consider the (generalised or weak) Laplace operator denoted by 



Hj(/2) -> H _1 (/2) 



i y s/u := — Au := — div Vu = — J2i=i 



u. 



(2.1) 



GRADIENT RECOVERY IN PARABOLIC PROBLEMS 



3 



We denote by L 2 (f2) the space of square summable functions on i?, with inner 
product and norm respectively defined by 

</,<?) := / /WsMdx and ||/|| := (fj) 1/2 . (2.2) 
Jn 

We will use the standard |Cia78| IEva981 ] Sobolev spaces 

H 1 ^) := {(j) G La(tf) : G L 2 (/2)} , (2.3) 

Hj(fl) i-I^GH 1 ^) : cf>\ dQ = Q} (2.4) 

and H _1 (/2) := dual(Hj(/2)) . (2.5) 

Let T > 0, the model parabolic problem consists in finding a function it G 
L 2 (0,T;Hj(/2)) and <9 t w G L 2 (0, T; H -1 (J?)) such that 

d t u(t) + ^u(t) = /(•,*), for all t G (0,T] , 

u(x,0) = ito(a;), for x G i?, (2.6) 

u(x, i) = 0, for (x, t) Edflx (0, T] . 

We consider the case where ito G L 2 (i?) and / G L 2 (0, T; L 2 (/2)) for which the 
problem (2.6| admits a unique solution [Eva98, ]. 

Problem ( |2.6[ ) is understood in the following weak form 

{d t u{t)A) + a{u(t),ct>) = {f{t),4>) V0GHj(f2), ie (0,T] 
u(-,0) = ito(-), 



where (•,•) is defined in (2.2) and a((f>,ip) :— V^). The form a(-, •) is clearly 
bounded and coercive, i.e., 

a(0,0 > a IMI? V^GHj(^), (2.8) 

where a = (1 + Cp) _1 and Cp is the Poincare constant. The bilinear form defines 
an inner product on Hq(.!?) and hence we can denote the energy norm ||-||^ := a(-, •). 

These observations justify our use of ||-|| a (instead of ||-|| H i ( fi J as the norm of 
Hg(j?) to be with the implied dual norm on H _1 (.f2) in (2.5). 

2.2. Spatial discretisation. Let 2? be a conforming, not necessarily quasiuni- 
form, triangulation of f2, i.e., (1) if G <^ means if is an open simplex (triangle for 
d = 2 or tetrahedron for d = 3) , (2) for any if , J G =5^ we have that if n J is a full 
subsimplex (i.e., it is either 0, a vertex, an edge, a face, or the whole of if and J) 
of both if and J and (3) Uatg-T K = f2. The shape regularity of ^ is defined as 
the number 

M(^):= inf ~ (2.9) 

where is the radius of the largest ball contained inside if and hx is the longest 
side of if. An indexed family of triangulations {2? n } n is called shape regular if 

fj, : = inf n.(3T n ) > 0. (2.10) 

n 

We will use henceforth the usual convention where h : fl — > R denotes the mesh-size 
function of 3^, i.e., 

h(x) := hg-(x) := max/i^-, and h n := ft,^™. (2-11) 

With a triangulation 2? as described above, and an integer p > 1 considered 
fixed in the sequel, we may consider the finite element space 

V := {<£ G Hj(/2) : <Z>| K G P p VA'G^}; (2.12) 
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and P k denotes the linear space of polynomials in d variables of degree no higher 
than a positive integer k. The spatially discrete finite element solution in V, is the 
function U : [0, T] V such that 

(dtU,$)+a(U,$) = (f,$) V#eV, 

(2.13) 

U(x,0) = U° :=n y u Q (x) Vxefl, K ' 



where 7T V : hi(Q) — > V is a suitable projector (or an interpolator if the data uq is 
in a higher regularity subspace of La(/2), e.g., £?-wise continuous) and (•, •) is the 
same as in (2.7) and (2.2). 

We will often write the scheme ( 2.13[ ) in its pointwise form 

dtU + AU = P f and 17(0) = U°, (2.14) 

where the finite dimensional space operator A : V — > V is the discrete Laplacian 
defined, through the Ricsz representation in V, by 

(AV,$) = a(V,$) V<£eV, (2.15) 

and P : L2(J7) — > V is the L2(^7)-projection operator such that, for each v € La(i?), 

(P v,$) = (v,$) V^eV. (2.16) 

The pointwise form is convenient as it allows for a more compact notation. 

2.3. Fully discrete scheme. Subdivide the time interval [0, T] into a partition of 
N consecutive adjacent subintervals whose endpoints are denoted to = < t\ < 
. . . < tN — T. The n-th timestep is defined as r„ := t n — t n —\. We will consistently 
use the shorthand F n (-) := F(-,t n ) for a generic time function F. 

The backward Euler method consists in finding a sequence of functions, U n € V™, 
such that for each n = 1, . . . , N we have: 

_ (u n -A n U n - 1 ,$) + a(U n ,$) = (/",<£) V^eV, 

T n (2.17) 

u° = n° UQ , 

where yl v : C°(f2) -> V denotes the La grange interpolation operator, A n := A , 
and 77° is defined as II y . 

Note our nonrestrictive use of the Lagrange interpolator as a "data-transfer" 
operator from a finite element space to the next. We do this to reflect exactly what 
we do in practical computations (where interpolation is faster than averaging) . All 
our analysis applies, however to a different data-transfer operator, including the 
h%{Q) projector, if necessary. 

As with the semidiscrete scheme the fully discrete scheme can be written in a 
pointwise form as follows: 

-+A Y U n = Pgf n and U° = 77°u , (2.18) 



where A n = A v " and P£ = P V " (cf. ( |2.15| ). 

2.4. Recovery a posteriori estimators. The stationary elliptic problem corre- 
sponding to a steady state of the evolution equation (2.6) is, 

given g € L 2 (i7), find w G Hq(J?) such that a/w — g, (2-19) 

where the operator is understood in a generalised sense and the solution is a weak 
one. The finite clement discretisation of the elliptic problem (2.19) consists in 

finding W G V such that a(W, $) = {g, <5) V <P G V. (2.20) 

We shall henceforth denote by w and W the solutions of (2.19) and ( 2.20[ ). 
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From the literature on elliptic a posteriori estimation [AO00, BR78l lCla781IVer96l 
BX03a, ZZ87, ] there is a variety of ways to compute upper and lower bounds for the 
error in some functional space X (e.g., Hj(f2), L2(J2) and L ocl (f2)). For instance, 
a generic upper a posteriori X -norm error bound takes the form 

\\w-W\\ x <£[W,g,X,V], (2.21) 

where S is an appropriate estimator functional. 

One way of providing an estimator functional consists, for example, in start- 
ing by applying a gradient postprocessing operator (postprocessor), say G, to the 
approximate solution W . And then proving that ||GW — VW\\ is equivalent to 
the error \\Vvu — VW\\. Gradient recovery operators form a subclass of gradient 
postprocessors. 

Recovery operators can be built in a variety of ways such as local weighted aver- 
aging (where the gradient is sampled from neighbouring elements) [PTc03, ], discrete 
L2(i?)-projection (using least squares fitting) [ZZ871 ] or global L2(i7)-projection 
(where a full discrete problem is solved) [BX03a, ]. In our numerical results we 



use local weighted averaging, defined explicitly in (5.4), but our theoretical results 
can be applied with any choice of recovery operator that provide upper and lower 
bounds for the elliptic problem. The fundamental idea behind these approaches is 
to build an approximation GW of Vro which is more regular than the piecewise 
discontinuous gradient VW; the extra regularity is aimed at obtaining a higher 
approximation order. 

2.5. Definition (gradient recovery operator, from |AO00j ). A gradient recovery 
(ZZ) operator on V is a linear operator G : V — > V d which enjoys the following 
properties: 

Consistency: we have, with A y : C°(J7) — > V denoting the Lagrange inter- 
polator, 

G(A y v)\ K = Vv\ K VuG P p+1 , K G ST. (2.22) 

Local bound: there exists a Czz > such that 

||GV|| LooW < Czz II Vf|| Loo (*) VUgV, Kef, (2.23) 

where K is the patch generated by K (the union of all L G & such that 
LC\K ^ 0). 

For simplicity, we assume that the operator is in a mesh-local relation with 
W noting, nonetheless, that global methods such as the global L2(i7)-projection 
proposed by [13x03a, BX03b] exist and can be included in our discussion. 

Under certain regularity assumptions recovery estimators are shown to be asymp- 
totically exact. For instance, |Zla77j shows that if w £ H S+1 (J?), with reference to 



(2.19) and (2.20), its approximation W € V satisfies the following superconvergence 
||V(W - A v w)\\ = 0(h p+<: ) for some C G (0, 1] . (2.24) 



property: 



A review of superconvergence results is given by |KN87| . If (2.24) is satisfied then 



the recovered gradient also satisfies the following superconvergence property [AO00, 
]: 

|| Vw - GW\\ = 0(h p+c ) for some Q G (0, 1] . (2.25) 
The reach of Zlamal's result is appreciated by stating the following consequence. 

2.6. Lemma (gradient recover y a p osteriori estimate from lAOOOj ). Let V be the 



finite element space defined in {2.12) and G : V — > V d a gradient recovery operator 
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according to \2.5 Ifw,W are the solutions of (2.19) and (2.20), respectively, and 



(2.25) holds then the recovery operator is asymptotically exact, in the sense that 

IIVW-GWII 

lim '' , jir = 1. (2.26 

hs-^o \\V(W-w)\\ K ' 

Thus, there exist S > 0, such that 6 (h) — > as h — > and 

(1 - S ) || W - GW\\ < \\V[W - w}\\ < (1 + So) ||VW - GW|| (2.27) 

for all partitions ST of ft satisfying h& < ho- 

2.7. Remark (recovery in absence of regularity). Lacking Zlamal's regularity as- 
sumption, recovery-based estimators are empirically observed to be efficient, reliable 
estimators, even on meshes with low shape-regularity |Car04[ ]. 

For more details about recovery-based estimators we refer to the available liter- 
ature |BX03al lBX03bl IXZfMl ILZMI bTOOOl IFVfiffl ]. 



2.8. Definition (gradient recovery a posteriori estimator functional). Lemma 2.6 
then justifies the use of the recovery estimator in the Hj(/2)-norm (and by equiva- 
lence the energy norm) by defining, for the rest of this paper, the gradient recovery 
a posteriori estimator functional 

&[V]:=&[V,tia(f2),V] :=||GV-W||, for V e V, (2.28) 

where G is a gradient recovery operator as defined in |2.5| 



2.9. Assumption (elliptic a posteriori error estimates). We will consider henceforth 
the blanket assumption that for a fixed ho, there are some Cq < Cq, such that 



for any V with mesh-size h < ho, for w and W solutions of (2.19) and (2.20), 
respectively and § defined in |2.8| the following bounds are true 

c S[W] < \\V[W -w]\\< C £[W] . (2.29) 

Optionally, we will assume asymptotic exactness, in which case 

C < l + B(h Q ) and c > l + (3(h ), (2.30) 

for some continuous functions B and j3 that vanish at 0. 



Assumptions (2.29) and (2.30) are true, modulo hierarchic oscillation terms of 
the data function g in (2.20). These assumptions are thus justified, for example, 
when g is in a finite dimensional space, for example g € V as we shall assume in 
the sequel, by isolating the bulk of the oscillations in data-approximation terms. 
For a thorough discussion of the oscillation in the context of recovery, we refer to 
[FVf)6] . 

The lower bound is not needed for the theory to be developed herein, as we will 
prove only upper bounds. Nonetheless, this property is required for the efficiency 
of the parabolic estimators in practical situations. 

3. Semidiscrete scheme 

To make the link between the parabolic problem and the elliptic recovered gra- 
dient estimates we utilise the elliptic reconstruction technique |.\1\03, [LM06 , ]. To 
make the discussion more accessible, we first do this for the spatially (semi) discrete 
scheme. We divide the error into two parts — one called elliptic error the other 
parabolic error — via the elliptic reconstruction of the discrete solution. Because 



the elliptic error can be directly bounded under the blanket Assumption 2.9 it is 
enough to show that the full error can be bounded in terms of the elliptic error 
only. This result is in accordance with the fact that the parabolic error on uni- 
form meshes is of higher /i-order in the energy norm with respect to the elliptic 
(and thus the full) error, as observed by [LM06 . The main result of this section is 
summarised in Theorem |3j3J 
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3.1. Definition (elliptic reconstruction). The elliptic reconstruction operator is 
denned as & : V -> Hjjj(/2) such that 

= (3.1) 



where A is the discrete elliptic operator defined in (2.15|. In weak form, equation 
(3.1) reads 

o(«7,#) = {AV,*) V#eH£(fl), (3.2) 
and it is well defined in virtue of the elliptic problem's well Sydney's. We will refer to 
the function 3%V as the elliptic reconstruction of V, while the elliptic reconstruction 
operator & will be called the reconstruction operator (or just the reconstructor) 
from V. 

If U(t) denotes the solution of (2.141 at time t, we shall indicate by uj(t) its 
reconstruction fflU(t). 

Thus, posing g(t) := AU(t), we then see U(t) is the finite element solution 
corresponding to the elliptic problem of finding u>(t) G H (J?) such that &tfu)(t) = g. 

3.2. The error and its splitting. For the whole of this section we shall consider 
u to be the solution of (2.6 1, understood in the weak sense, and U its semidiscrete 
approximation given by (2.14). The corresponding semidiscrete error is defined by 

e(t) := U(t) -u(t), (3.3) 



and can be split, using the elliptic reconstruction uj = MU defined in §3.1| as 
follows: 

e(t) = (w(t) - u(t)) - (w(t) - U(t)) =: p{t) - e(t). (3.4) 

We shall refer to e and p here defined as the elliptic (reconstruction) error and the 
parabolic error respectively. Using this notation we have the estimate 

||V[tf-u](t)|| < ||Vp(t)|| + ||Ve(t)||, (3.5) 



where, following the remarks made in Dcfinition |3.1| and Assumption |2.9[ the elliptic 
error can be bounded by the computable elliptic a posteriori estimator functional 

S: 

l|e(*)||„ = ||Ve(*)||<Co^(t)]. (3.6) 
It is therefore sufficient to bound the error's energy norm using the elliptic error's 
energy norm. 

3.3. Lemma (elliptic energy bound for parabolic semidiscrete error). If e, e are 
defined as in |3.2| then, for each is [0, T], we have 



Ht)\\ + l|e(*)||;da<||e(0)ir+/ \\e(sX + 2 (PJ(s) - f(s),e(s)) ds. (3.7) 
Jo Jo 

Proof From the the exact problem (2.6), the semidiscrete scheme ( |2.14[ ), and the 
splitting (Q 

d t e + ss/p = d t [U -u] + s^[uj - u] = d t U + AU - d t u -srfu= P f - f. (3.8) 
Testing with e we obtain 

(d t e,e)+a(p,e) = (P f-f,e) (3.9) 

and thus 

^d t || e || 2 + ||e||^ = (F /-/, e )-a(e, e ). (3.10) 
Integration from to t yields 

l|e(t)f+2 / \\ef a = ||e(0)|| 2 +2 f (P f - f, e)-2 f a(e, e) V< £ [0, T] . (3.11) 
Jo Jo Jo 
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Hence, by Young's inequality on a(e, e), we have 



l|e(*)ir + 2/ ||e||^<||e(0)|r + 2/ (PJ /, e) + ||e|£ + / HC> (3.12) 
Jo Jo Jo Jo 

whereby the claim is verified. □ 

3.4. Remark (proliferation of \[2 syndrome). Let a, b, c > such that a 2 < c 2 + ab, 
then, by Young's inequality, it follows that a 2 < 2c 2 + b 2 . Note however that the 
factor "2" in 2c 2 is not needed, in that we also have that a < c + b. If 1 > a ~ 
c ^> b > 0, then the first bound provides a/c « -\/2 whereas the second bound gives 
a/c « 1, which is tighter. 

The following result, which generalises a < c + b, is extremely simple yet useful 
in avoiding this "proliferation of \[2 syndrome" from repeated usage of Young's 
inequality. 

3.5. Proposition (L 2 simplification rule). If a, b £ R , N E N, c € R and f,g € 

L2(D), for some measurable domain D, are such that 



|a| 2 + ||/|| 2 <c 2 +aTb+ / /«?, 

JZ5 



then 



l/ll 



1/2 



< 



|c|+(|t 



1/2 



(3.13) 
(3.14) 



where all the vector norms are Euclidean, and the function norms L/2(-D). 



Proof 

Denote by a := (|a| , ||/||) and [3 := (|b| , ||. 9 ||). 

If \at\ < \(3\ then (3.14) is trivially satisfied. Otherwise we have \a\ > \f3\ whereby 
(3.13) and the Cauchy-Bunyakovskiy-Schwarz inequality imply that 



lal 2 <c 2 



< c 2 



lallbl 



|/|||MI + |i9|(|a|-|/3|) 



2\ a \\/3\-\(3\< 



Hence (|a| — |/3|) 2 < c 2 , and thereby 



lal < c 



1/31, 



as claimed. 



(3.15) 



(3.16) 
□ 



3.6. Theorem (a posteriori semidiscrete error estimate). With u and U as defined 
by {2.6) and [2.13), respectively, and an estimator functional S as defined in [2.28), 
we have 

\\u(t)-u(t)f + y o iie/-<J 1/2 

< \\U(0) - u(0)\\ + Co IW]|| L2[0 ,t] + 2 \\Pof - /IIl 2 (o,t ; h-m«)) ■ (3-17) 
Proof Using Lemma [3~3| we have 

(3.18) 



e(i)|| 2 + / || e || 2 <|| e (0)|| 2 + / N| 2 +2/ (Plf-f,e) 



Using Proposition 3.5 we obtain 



|e(t)|| 



1/2 / r t xl/2 / f t 

<(j|e(0)|| 2 + J o Ml) +2^jf \\P f-ff^ 



1/2 



(3.19) 



Assumption (2.29) and the discussion in ^3.2 ensure then that 
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ft \l/2 

\e(t)f+ / ||e|| 2 



r t \l/2 / r t \l/2 

|2 , r~a I am? \ : r, ll r> £ j- 1 1 2 



<(||e(Q)|r + £7^y +2M ||P /-/||^-i (fl )] , (3.20) 

which implies the claim. □ 

3.7. Remark (short versus long integration times). The bound for the pointwise 
in time L 2 (.!?) error, ||e(t)||, appearing on the left-hand side of (3.17), is tight only 
for very short times. For example, it is well-known that on a uniform mesh of size 

a-t 2\ 1 / 2 
||e|| a ) is 0(h p ), while ||e(i)|| is 

0(^ +1 ). 

3.8. Remark (dealing with the H _1 (i7) norm). In practise the H _1 (X?) norm can 
be well approximated as shown by Lemma[3.9| so, in the lack of a priori information, 



the last term in (3.17) may be replaced using the Poincare inequality 

2 Iko/ _ /IIl 2 (o,t ; h- 1 (^)) - 2( ^p(^) \\ p of _ f\\h 2 (nx(o,T)) • ( 3 - 21 ) 

It is also possible to obtain bounds by using the Cauchy-Bunyakovskiy-Schwarz 
inequality for L 2 (J?) on the term (P f — f, e) — rather than the (H _1 , Hj) duality — 
and "absorb" the resulting sup[ t ] ||e|| into the first term on the right hand side of 
(3.17). However, whenever possible, we shy away from this procedure as it incurs 
in artificially higher constants and a Li [0, T] accumulation on the right-hand side 
while the energy term on the left-hand side accumulates like L 2 [0, T]. This time- 
accumulation disparity between the error and the estimator is likely to result in 
an error-estimator ratio bound that has the order of vT, that is, although having 
the right order of convergence, the estimator will overestimate the error over long 
integration times. 

We show now how to practically approximate the H _1 (i7) norm of an arbitrary 
given function v G L 2 (J?). 

3.9. Lemma (computing the H" 1 (f2) norm). Let v € L 2 (,!?), consider the functions 
i> e h£(J7) and & E V such that 

*/ip = v and A<P = P a v, (3.22) 

where A and P are the discrete Laplacian and the L 2 (/2) projection on V, respec- 
tively. Then, recalling our convention whereby IMIh^j?) = \\^ v \\ we have 

NlH-Hn) = l^llHj(fl) = II V- - nl hi o) + Wnlun) ■ (3- 2 3) 



Proof With if) and & as given in (3.22 1 we have $eV 

{sitil)-AV\$) = {v-P fJ v,$)=Q l (3.24) 
i.e., that ip — $ is Galerkin-orthogonal to V. Also, we have 

IMIh-Vu) = Uhh(n) ■ (3- 2 5) 

Indeed, on the one hand 



sup ^ ^ = sup V ^ 

M^ n) II^IIhiw _ 



(3.26) 
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and, on the other hand 

By the above, Galerkin-orthogonality and Pythagoras's Theorem, we have 

IMIh-i(«) = II^IIhjc^ = - + ll*fcn) • (3-28) 

□ 

3.10. Remark (the H~ (J?) norm approximation). The next-to-last term \\ip — ^|| H i(™ 
is the error of a function and its Ritz projection. This can be easily estimated with 
a fully computable a posteriori estimator functional S such that 

||^-*llHS(fl)<*[*.».V] = 0(/i5), (3.29) 

where hy is the "mesh-size" of the space V. 

Hence the H _1 (r2) can be computed using the relation: 

\Hl-Ho) = + C[#, «] 2 , where ([9, v] < g[W\ . (3.30) 

The term ||^|| H i(fi) is clearly computable, by computing W, which involves one 
L 2 (i7)-projection, one stiffness matrix inversion and one (discrete) energy norm 
computation. Furthermore 

H«llH-(fl) = ll^(fl) + 0(^)- ( 3 - 31 ) 
Hence, if \P is finite with respect to the mesh-size hy, i.e., H^l = O(hy), then we can 
approximate the H _1 (i7) of a function with as much precision as the finite element 
method allows it for the energy norm. On the other hand if \P is small, precisely, 
W = 0(hy) with s > (implying that v is small as well), then this result has to 
be handled with more care for the error to be of some order of hy higher than the 
computed quantity. 



3.11. Remark (sharper versions of Theorem 3.6). The error estimate (3.17) can be 
tightened further to 



^l|e(t)|| 2 + / ||e 112 



1/2 



< ^ ||e(0)|| + (J* \\P Q f - f\\^ {Q) + Clspf^ 2 - (3-32) 



But this estimate becomes noticeably better only when one of the terms ||e(0)|| or 
ll-Po/ — fWa-^fa) dominates the &[U] term, which should not be allowed to happen. 
So there is no need to lengthen the discussion by insisting on such tight bounds, as 
long as it is possible to obtain the elliptic a posteriori estimate constant Cq in the 
leading term on the right-hand side. 

4. Fully Discrete scheme 
The main result of this section — and the paper — is the a posteriori error bound, 



stated in Theorem 4.6 on the error between the approximate solution U of the fully 



discrete problem (2.18) and that of the exact problem (2.6). 

The analysis in this section follows narrowly the one we performed in £j3j albeit 
with the complications that the fully discrete scheme imports. We will first extend 
the discrete solution sequence to a continuous-time function. Then we derive an 
error identity on which we mimic the energy techniques of fj3]to bound the error's 
energy norm in terms of some residual terms and the elliptic error's energy norm, 
which is finally controlled via gradient recovery estimators. 



GRADIENT RECOVERY IN PARABOLIC PROBLEMS 



11 



4.1. Time extension of the discrete solution. Recalling the fully discrete 
scheme (2.18), the fully discrete solution is the sequence of finite element func- 
tions U n £ V™ defined at each discrete time t n , n = 0, . . . , N. Define the piecewise 
linear (affine) extension 

N 

U(t):=J2U n ln(t), (4.1) 

71=0 

where we use the one-dimensional piecewise linear continuous Lagrange basis func- 
tions , defined for t > , as 

!{t-t n -i) /r n , for t £ (t n -i,*n] (and n > 0) , 
(t„+i - t) /r n+1 , for t £ (t n ,t n+1 ] (4.2) 
0, otherwise. 

We warn the reader that we use the same symbol, U, to indicate the fully discrete 
solution's extension to [0, T], as the one we used for its semidiscrete counterpart in 

§s 

4.2. Elliptic reconstruction and error splitting. Next we define the elliptic re- 
construction, needed for the following analysis, similarly to that of the semidiscrete 
sche me (cf. (3.1 )). For each n £ [0 : N], with the discrete elliptic operator A n as in 
2.3 we define the corresponding elliptic reconstruction operator & n : V™ — > Hq(^2), 



for each V £ V™, by solving for 3% n V the elliptic problem 

i/J n y = A n V, (4.3) 
which can be read in the weak form as 

a{M n V : <P) = (A n V, <P) V#eHj(tf). (4.4) 

We denote 

oj n := M n U n , for each n = 0, . . . , N, (4.5) 

and this sequence's piecewise linear extension by u : [0, T] — > H (J2), i.e., 

N 

u{t):=J2" n ln(t)- (4-6) 

As in the semidiscrete analysis we introduce symbols for the full error e := U — u, 
the elliptic error e := to — U and the parabolic error p :— uj — u, whereby 

e = p-e, (4.7) 

and, based on the Assumption |2.9| 

IKi)L < c £[u n i n (t) + u n -H n ^(t)} 

<C (^[C/ n ]Z n (f) + ^[[/ n - 1 ]/ n _ 1 (f)) forte [t n -i,t n ]. 

The last step is guaranteed by the linearity of the operators G and V, hence the 
homogeneity <f[AV^] = |A| \\GV — W||, and by the triangle inequality . 

4.3. Lemma (parabolic error identity). For each n — 1,...,N and each t £ 
(t n -i,t n ) we have 

d t e(t) + *p(t) = (A n U n ~ 1 - C/"- 1 ) /r„ + sf[u{t) - ut n ] + P n / n - f(t). (4.9) 

Proof By the definition of U, ( |4.1[ ), for each n = 1,,..,N and t £ (t n -i,t n ) we 
have 

d t U(t) = U n l' n (t) + U n - 1 l' n _ 1 (t)=(U n - t/"- 1 ) /r„ (4.10) 
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and using the fully discrete scheme (2.18), we have 

dtU(t) + = {A n U n ~ l - U"- 1 ) /r n +(U n - A 71 !/ 11 - 1 ) /r„ + A n U n 



(4.11) 



= (A n U n - 1 - U 11 - 1 ) /r„ + P "/ n . 
Hence 

d t U{t) + s^uj{t) = (A 11 ^- 1 - U n ~ x ) /r„ + j*[w(t) - uj n ] + P 7l f n (4.12) 
and, using the exact PDE (2.6), we get 



d t e(t) + sfp(t) = d t U(t) + gfu(t) - d t u{t) - srfu{t) 

= {A n U n - x - C/"- 1 ) /t„ + tf[u{t) - u n ] + Pgr - f(t), 
as stated. □ 



(4.13) 



4.4. Definition (a posteriori error indicators). The notation we introduce here will 
be valid for the rest of the article. 

elliptic error indicator via recovery: 

e n := £[U n ^l{C2),V n ] = C \\VU n - G n [U n }\\ , (4.14) 

with the functional S as defined in fJIU ancQ 

4 := \{4 + 4-1 + e»£n-i) < 2& + e «-i) ' ( 4 ' 15 ) 
time-discretisation error indicators: 

1 / \\Pof n ~ A n dU n - (P " _1 /" _1 ~ A n ~ 1 dU n ~ x ) \\ a -i, n) for n > 2, 
6n ^ ^ lll^o 1 / 1 -^at/ 1 - A (7°|| H _ 1(r2) for n = l, 

(4.16) 



where (?{/" := (Z7 n — f/ n : ) /t„, (cf. Lemma 3.9), also possible to use in 
its alternative (faster to compute but not as sharp) version 

e n :=C^\\U n - 1 -U n \\ a , (4.17) 

where C M is dependent on the shape regularity ji of the family of triangu- 
lations defined in (2.10). 
mesh-change (coarsening) error indicators: a main mesh-change indi- 
cator 

7n := r- 1 \\A n U n - x - U n - X \\^ m , (4.18) 
and a higher order mesh-change indicator 

- _ r l j\\h n (P^f n -A n dU n -P^- 1 f n - 1 +A n - 1 dU n - 1 )\\, n>2, 

7 " llM^-^W-^t/ )!!, n = l, 1 ' 

where h n (x) — max{h n -i(x),h n (x)} for x £ Q and a constant C'. 
data approximation error indicator: 

Pn ■■= r- 1 I" \\PSf n - /(t)|| H -i(fl) d*. (4.20) 

Jtn-l 

4.5. Remark (computing H _1 (J7) norms). Clearly the H _1 (]7) norms appearing 
in Definition [4T4] cannot be computed in practise. The corresponding indicators can 
be replaced by upper bounds using the (dual) Poincare inequality 

Mh- 1( o) < Cp WW ■ (4-21) 
Other alternatives, as described in Remark 3.8 are possible but will not be described 
here. 



4n the numerical experiments we use (e^ + e'f instead of s n - 
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4.6. Theorem (a posteriori estimate for fully discrete scheme). Let the sequence 
(U n ) ne[0:N] , U n £ V", be the solution of the fully discrete problem ( |2.17[ ) and U its 



problem (2.6) then 



piecewise linear time-extension as in (4.1 1. Let u be the exact solution of the exact 

\l/2 



\U N -u(T)\\ 



+ / \\U(t)-u(t)\\ a dt\ < 



\U(0)-u(0)\\ 
V2 



Vn 



(4.22) 



where the (global) error estimator is given by the following discrete L2(0,T) sum- 
mation of the error indicators defined in 94. 4\ 



N 



Vn 



(4.23) 



Proof The proof shadows that of Lemma |3.3| and Theorem |3.6[ but we must take 
into account the complications arising from the time discretisation. For the reader's 
convenience we divide it into steps. 



Step 1. Using the notation from Lemma 4.3 and identity (4.9 1 therein we have 
that 



d t e(t) + £/e(t) = ^e{t)+{A n U n - 1 - U^ 1 ) / 



+*/[u>(t) -cJ n ) + P "/ n -/(*)■ 



Testing this with e we obtain 



1 



^ d t ||e(t)|| 2 + \\e(t)\\l = a(e(t),e(t)) + ((A^ 1 - C/"" 1 ) /r n , e(t)> 
+ (^[«(t) - u> n ] , e(t)> + {PSf n ~ /(*), e(t)> , 
for all t € (tn-i, t n ) and each n = 1, . . . , N. Integrating over [0, T] gives 



JVI 



/2+ / ||e(t)||^dt = ||e°|| 2 /2+ / a(e(t),e(t)) dt 



o 

N -U 



/ ((^^-tP-^/rn.eft)) 



+ a(u(f) - w n , e(t)) + (P "/ n - /(*), e(t)) dt 



(4.24) 



(4.25) 



(4.26) 



0! + S 2 + ^ 3 + ^4 + ||e°|| /2. 
We proceed by bounding each of the terms Bj, j = 1,...,4, appearing in the 



right-hand side of (4.26) 



Step 2. The first term to be bounded in (4.26) yields the spatial discretisation 
error indicator as follows: 



pT N pt n 

Bi= / a(e(t),e(t)) dt = J2 a(e(t),e(t)) dt 

Jo n=1 Jt n -i 



N 



> I / IN 

n=l \ Jt n-l 



1/2/ pt 



.1/2 N 
2 \ „ ~ _l/2 



< 



n=l 



tn-l 



1/2 



(4.27) 



where we have used (4.15) and in view of (4.8) and (4.15), we may write 



e|| Q < E n-i / ^n-i + 2e n _ie„ / l n -il n + e 



I 2 - e 2 r 



(4.28) 
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The second term in ( |4.26 ) contains mesh-change term which we bound as follows: 

N f-tn 

B * = E / ((^ n U n - 1 - U 11 - 1 ) /r n ,e(t)) dt 

n=l Jtn-l 



N 



1/2 



(4.29) 



tn 



n=l 
N 

<E^« /2 

n=l 

where 7„ is defined by ( |4.18 ). 

Similarly the data error term is bounded as follows 

i-T N / /•«„ 



64=/ (P n r-/W, e W>d<<f]/3„ry 2 f /" \\e\\l 
Jo n=1 \it„_i 



1/2 



(4.30) 



where f3 n is defined in (4.20). 



Step 3. The third term in (4.26 1 yields a time discretisation term and is a bit more 
involved to estimate. Using the definition of uj 7 \ w and 3% n ', given in (4.3) and (4.6), 
we observe that 

N f-tn 

63 = E / a(u-oj n ,e(t)) dt 

N rU 



N 



(4.31) 



.1/2 



n=l 



t„ \l/2/ /■*•" 

2 ' 



1/2 



= V / a^^M 11 - 1 ^- 1 + l n {t)M n U n -^ n U n ,e{t)) dt 
= J2 / In-i^a^ 1 - 1 ^ 1 - 1 -M n U n ,e{t)) dt 
= J2 I ln~i{t) {A n ~ l U n ~ l - A n U n , e(t)) dt 

I 2 

l n-l 

v ' V^t„_i 

E^jT ||e|| 2 a 

„=1 V 

where in the last passage we use the discrete scheme (2.18) for the substitution 

A n U n = (A n U n - 1 - U n ) /t„ + Pof n for n > 1. (4.32) 
Step 4. Grouping together ( |4.26[ ), ( |4.27[ ), ( |4.29[ ), ( |4.30| and ( |4.31[ ), we obtain 

^ |e(t)||2 dt 



Kir/2- 



^ ll e °ll / 2 + E^« + ^ i+/3 « + 6, ") r ' 

n=l 

Using an L2 simplification (cf. |3.5[ ), we conclude that 

|2 „T \V a 



1/2 



1/2 



/|| e JV||2 „T \V II ON /W \l/2 

Li + ^ ||e(t)||* dij < iLJi + I J2(e n + 7« + /3„ + 0„) 2 rj . (4.34) 



(4.33) 



Referring to the notation in (4.1 ) and Definition 4.4 we obtain the result. 



□ 
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4.7. Remark (the alternative time indicator). Assuming there is no mesh change 
from time t n -\ to time t n , then the discrete Laplacians A 11 ^ 1 and A n , defined 
in (2.15), coincide. Thus the time discretisation error indicator 9 ni which is part of 



the estimator rjpj in Theorem |4.6[ can be written as 

^ = ^ W AnUn - An ~ lun ~ l L-Ho) = % W And ^W^{o) 



(4.35) 



In the form given in (4.16) and using the dual Poincare inequality (4.21), this 
indicator is easily bounded. 

A more precise, but slightly more expensive, calculation can be done using 
Lemma 3.9 The same idea, will be used in the next result where we show that the 
indicator 9 n is equ ivalen t, up to higher order terms, to the alternative time indica- 
tor 9 n , defined in (4.17), which requires only an energy norm computation. This 
alternative time indicator, which is more common in energy estimates [Pic98[ e.g.], 
9 n is also more "natural" , as it measures the time derivative in the energy norm 
as opposed to the H _1 norm of the time derivative of AU. Due to mesh-change 
effects, this simpler indicator comes at the (affordable) price of having to add the 
higher order mesh change term 7„ to the otherwise simpler 7„. 

4.8. Theorem (alternative time estimator). With the same assumptions and no- 
tation of Theorem 14.61 we have 

\l/2 



U 



N 



u(T) 



\\U(t)-u(t)\\ 2 a dt 



< 



\U(0)-u(0)\\ 
V2 



(4.36) 



where the (alternative global) error estimator is given by the following discrete 



L 2 (0, T) summation of the error indicators dehned in { 4.4 

- In + Pn + 9n 



N 



n=l 



(4.37) 



Proof We proceed similarly to the proof of Theorem |4.6| in steps. The notation is 
the same and steps 1 and 2 are identical. 

Step 3. This step starts similarly to its homologue in the proof of Theorem |4.6| by 
observing that 

N 

l n -x{t) {A n - l U n - 1 -A n U n ,e{t)) dt. (4.38) 



n=l 

The function A^U^ 1 



A n U n belongs to 



but in general it is in 



neither of V™ nor V n " 
Clement-Scott-Zhang interpolator denoted respectively by 



Thus, to proceed, we use the L 2 (^7)-projection and the 



P n : L 2 (f2) 



and 77™ : L 2 (J2) -> V" HV r ' 



(4.39) 



We recall that the operators 7J n and P n are both known |SZ90( ICar02[ resp.] to 
enjoy the following stability properties for all n = 0, . . . , TV: 



n n (, 
P n 6 



< | 



V^eH 1 ^), 
V^eH x (/2), 



(4.40) 
(4.41) 



where fi is the shape-regularity of the triangulation family {^ n } n=0 N defined in 
(2.10). Furthermore, the following interpolation inequality is valid |LM06( §B.3] 



(^-77"V) fh n \\ <C^U\\ a WeHj(^),n = l,...,TV, 
where h n :— max {h n , /i„_i}. 



(4.42) 
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Step 4. Using these operators, we derive that 
N r t„ 

B 3 =^2 {A n - x U n ~ l - A n U n ,P n e(t))l n - 1 {t)dt 



= V f n ((A n - 1 U n - 1 ~A n U n ,P n e{t)-fl n P n e{t)) 

+ (A™- 1 !/* 1 - 1 - A n U n ,n n P n e{t))y n - 1 {t) At 

<V / (\\h n (A n - 1 U n - 1 -A n U n )\\\\h- 1 (P n e(t)-n n P n e(t))\ 

+ a(C/"- 1 - U n , fl n P n e{t))Ya-i{t) dt. 
Using inequalities (4.401, (4.41) and (4.421, we get the bound 

N f-tn 

B 3 < T / (C 3)M fkniA"- 1 ^- 1 - A n U n )\\ \\P n e(t)\\ a 

n=1 J tn-1 



(4.43) 



<Y,( c ^f h ^ An ~ lu 

n=l 



C ltll DC/"- 1 - U n \\ a \\P n e(t)\\ a )l n ^(t)dt 

n-lrjn-l _ ^[/™)|| + £ Nj/n-l _ [/« I 



(4.44) 



X G 



241 



|e(t)|| in-i(t)dt 



t„_i 



< 



by taking := C Wl C 2 , M /3, C/ := C^C^JZ in and gJ9) for the last 

step. 

We may now conclude exactly like the last step in the proof of Theorem 4.6 albeit 
with n replaced by j n +d n . □ 

5. Computer experiments: convergence rates 

In this section and in £[6] we study the numerical behaviour of the error indicators 
and estimators and compare this behaviour with the true error on three model prob- 
lems. The C code that we used includes the adaptive FEM library ALBERTA [SS05, 
]. The quadrature formal error is made negligible with respect to other error by 
using overkill quadrature formulas (exact on polynomials of degree 17 and less). 

5.1. Benchmark problems. Consider three benchmark problems, the solution of 
which is known. Namely, take d = 2, each problem's data /, uq is then chosen such 
that the exact solution to |2.6| is given by: 

(5.1) 
(5.2) 



(5.3) 



i(x, t) = sin(7rf) exp^— 10 |a;| 2 ^ , 
i(x, t) = sin(20vrf) expf-10 \x\ 



, , . 2 arctanfe/xi) 0/5 / —1 

u(x. t) = t sin v ' J x 2/3 exp n 

y ' 3 \l-\x\ 2 



The domain Q for Problems (5.1) and (5.2) is the square S := (—1,1) x (—1,1). 
Problem (5.3), whose solution's gradient is singular at the origin, is considered on 
the L shaped domain Q = S x [0, 1] x [—1,0]. The benchmark problems (5.1) and 
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(5.2) have been chosen such that they can be compared with previous numerical 



studies LM06, 



For all Problems ( 5.1 H 5.3 1, we take zero initial condition, u — to avoid 



dealing with the initial adaptivity which is a side issue here. 



The solution (5.2 1 has a time dominant discretisation error, while (5.3) was 
constructed to have a dominant spatial error. It is the product of a linear function in 
time, a well known solution to Laplace's equation producing the spatial singularity 
and a mollifier. 



Problem (5.1) is used to test asymptotic behaviour of the indicators under uni- 



form space-time refinements further in S5.5 Problems (5.3) and (5.2) will be used 
to test the adaptive strategies in Sj6j 

5.2. Gradient recovery implementation. The recovery operator, G n , is ob- 
tained by taking the discontinuous gradients of the numerical solution at the super 
convergent sampling points [AOOOl ] (and references therein). The recovery op- 
erator used here is built in the following way: fixing V E V n , for each degree of 
freedom x, we define 

G"[V](x) := W ^\ K[X) > (5-4) 

This defines a unique piecewise polynomial field G n [V] € V d . (Note that for- 



mula ( 5.4 ) is non trivial for only for those DOF that are are on the boundary of an 
element; for the internal DOF, that arise in using P p elements for p > 3, it is not 
necessary to calculate anything.) 

5.3. Definition (experimental order of convergence). Given two sequences a(i) and 
h(i) \ 0, i = I, . . . „ we define experimental order of convergence (EOC) to be the 
local slope of the loga(i) vs. logh(i) curve, i.e., 

mr/ , -s log(a(i + l)/a(i)) , , 

E(JL(a, h :i) ;= - — t— - w , ■ (5.5) 

log(fc(t + l)/fc(*)) 

5.4. Definition (effectivity index). The main tool deciding the quality of an esti- 
mator is the effectivity index (EI) which is the ratio of the error and the estimator, 
i.e., 

EI(t n ) := T) n /\\U- u|| La(0itB . H i (n)) . (5.6) 
If EI(t„) —¥ 1 as sup.,, n h n (x) — > then we say the estimator is asymptotically exact. 

5.5. Indicator's numerical asymptotic behaviour. In the following conver- 
gence rate tests we discuss the practical realisation of Theorems |4.6| and |4.8[ to 
which we refer for notation. 

We use a uniform timestep and uniform meshes that are fixed with respect to 
time. Hence for each test we have V™ = V° = V and t„ = r(h) for all n = 1, . . . , N. 
For each test we fix the polynomial degree p and two parameters k, c and then 
compute a sequence of solutions with h = h(i) = 2~ 1 / 2 , and r = ch k for a sequence 
of refinement levels i = I, . . . , L. 

Due to the finite element space invariance in time, the coarsening indicator j n 
vanishes and is thus not computed (this indicator will be discussed in fj6j. 

The initial value being zero makes the initial error U(0) — u(0) zero. Thus we 
do not need to calculate this term in the estimator. 

For all solutions the boundary values are not exactly zero, but of a negligible 
value, hence little interpolation error is committed here (nonetheless some care is 
taken when dealing with very small errors). Finally, the data approximation error 
term, j3 n , though important for highly oscillatory data, will not be studied here 
given the regularity of our data. 
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Therefore, what we compute on a space-time uniform mesh are the indicators 



e n and 6 n (or #„,7„), defined in [4.4 and the corresponding cumulative indicators 



B n)n=i,...,iV and ( 6) ™)n=i,...,Af defined by: 

m 

E fn ■ = ( £ » + £ n- 1 ) T « / 2 ( f ° r S P aCe ) ' 



n=l 



m _ (5-7) 

and e 2 rn := ^ 6£r„ or ^ (fi£ + ^ Tn (f or time) . 



n—1 n—1 



From the Theorems 4.6 and 4.8 we know that 



- u{t n )\\ 2 < El + 2 m + fcr n . (5.E 

71=1 

Our results and the comments are reported in the captions of figures. 



In Figures [1[[4| we visualise the results and comment them, for Problem (5.1) 
for conforming finite elements of polynomial degree p = 1,...,4, respectively. 
Having fixed p,k,c such that r = ch k , for each level i, we plot m and E m , 
\\U — u|| L2 ( 0t -h 1 ^))' their experimental order of convergence and the effectivity 
index FI(t m ) versus (discrete) time t m = 0, ...,T. The conclusion is that the esti- 
mator is sharp and reliable, but to achieve asymptotic exactness (or close) the time 
indicator must be made smaller than the space indicator by taking r <C h p . In all 



these tests we used the first form for O m appearing in (5.7). 

In Figure [5] we summarise a comparison between the two time indicators 9 n and 
ni showing that the latter yields a much sharper bound, but with the added cost 
of having to compute the higher order term j n . 

6. Computer experiments: adaptive schemes 

We present now an adaptive algorithm based on the error indicators defined 
in j|4.4| A s with many adaptive methods for time-dependent problems |Pic981 
SS05, CJ04 ( ], we perform space and time adaptivity separately. Adaptivity is 



controlled via the indicators rj n and rj n (or rj n ) — see Theorems 4.6 and 4.8 — which 
are kept under a given tolerance tol. 

Namely, at each timestep t n -i — ► t n , we use adaptive schemes for elliptic prob- 
lems as to minimise the indicators s n and j3 n . There are different strategies to 
perform the timestep adaptivity, all geared towards minimising 9 n (or var#„). Fi- 
nally, the coarsening estimator j n is minimised by precomputing it and performing 
only one coarsening operation at the beginning of each timestep. 

Note that it is not in the scope of this paper to prove any rigorous result about the 
adaptive algorithm and, based on heuristic arguments only, we use it for illustration 
purposes. 

6.1. Space adaptivity via maximum strategy. At each timestep an elliptic 
problem is solved. For linear elliptic problems, convergence of adaptive schemes 
is reasonably well understood IMNS02I [BDD04I ] so we follow the criteria given 
therein, namely the Maximum Strategy. 

The algorithm we used can be pseudocoded as follows. 

6.2. Space Adapt. 

Require: (f/ old , V old , tol £ , fc max , t, r, £, tol 7 ) 



Ensure: (U new ,V aew ) solution of fl2T7 ) 
procedure Coarsening 



7 = {l K )Kasr ■= Coarsening Preindicator(J7old, V old ) (cf. d A.12[ ). 
£T ■= Mesh(V old ) 
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find <g C ST such that J2k ^ il K f < tol 7 
ST := Coarsen (^."T) using fSS051 §1.1.2-1.1.3] 
end procedure 

procedure Maximum Strategy Refinement [SS05 
k := 



compute e n using (4.14| 



& := > refinement set 

while e n > tol e and k < fc max do 
for K G ,T l do 

if eL„ > i max Le ,y-. e| „ then 

^ := {if} U ^ > mark K for refinement 

end if 
end for 

ST := Refine(^,^) using [SS05 , §1.1.1] > hence upda te ([/ old ,V) 



set yTl/"- 1 := C/ old , r n = r, t„ = t and solve for U n in ( |2.18 l 
[/ := U n 



compute e n using (4.14) 
k := k + 1 
end while 
end procedure 

return (U, V) 

6.3. Coarsening. In time-dependent problems mesh coarsening , which is not to 
be confused with the coarsening needed in proving optimal complexity for adaptive 
schemes [BDD04, ], is used to reduce DOF that become redundant in time. 

Mesh coarsening is a delicate procedure and should be used sparingly as to avoid 



needless overhead computing time. In Algorithm 6.2 coarsening is performed only 
once, at the beginning, for each time-step. 

The coarsening strategy we propose is based on predicting the effect of a possible 
removal of degrees of freedom. The reason for this is that in ALBERTA (and many 
other finite element codes) upon coarsening, all DOF-dependent vectors (encoding 
finite element function coefficients) are "coarsened" via interpolation. This makes 
it possible to compute the effect of coarsening, and the coarsening estimator 7„ 



defined in (4.18), before mesh-change occurs. The details of this procedure are 
discussed in § |A 

6.4. Timestep control. Timestep control can be achieved using two different 
strategies. 

An implicit timestep control strategy used is ready implemented in ALBERTA 



SS05, ] using Algorithm 6.2 upon each timestep. 

Here we propose an explicit timestep control strategy which we have implemented 
in ALBERTA. The reason for this is that the implicit strategy, though better in terms 
of timestep determination, is very time-consuming as it requires the repeated solu- 
tion of the timestep. In contrast, the explicit strategy has a rougher — nonetheless 
still satisfactory — control over the timestep, but it is much faster. The conclusion 
is that the ideal control should be a smart implicit/explicit-switching algorithm. 

The explicit strategy can be described as follows. 

6.5. Explicit Timestep Adapt. 

Require: (r , t , T, J°,u°, tol e , & max , £, toly, t ol 9 , min , tol e ) 



Ensure: (r„, V™, U n ) n =x,...,N satisfying (2.17) and possibly / Q T \\U — u\\ 2 < tol 2 



(U°,y°) = Initial Space Adapt(^°, u°, fc max , £, k) > data interpolation 

n := 1 
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f~n • Tfi—\ 

while t n < T do 

(IP, V») := Space Adapt(f/«-\ V"" 1 , tol e , fc max , r B> f B> £, tol 7 ) 
compute 9 n 
if 6>„ > tolg then 
T~n+1 ■= T n j v2 
else if 0„ < tole,min then 

T n +1 ■— V^Tn 

end if 

tn+1 '■— tn + T n +1 

n := n + 1 
end while 

return [U n )n=x,...,N, 
where the global tolerance tol is given by the relation 



tol 2 = T(tole+tolg+tol 2 



(6.1) 



Note that this algorithm does not guarantee reaching a tolerance, unlike more 
sophisticated ones found in the literature |CJ04[ e.g.], but it guarantees termination 
in reasonable CPU times. 

6.6. Numerical results. In Tables [l}[3] we compare the implicit timestep control 



strategy described by algorithm 6.5 with a uniform timestep scheme. For the uni- 
form strategy we take a stationary mesh in time and set r = 0.04h 2 . We calculate 
the error for various numerical simulations using differing values of h using the 
uniform strategy and set those values as tolerances for the adaptive scheme varying 
£ appropriately. 

Each column displays results for either the uniform strategy or the adaptive 
strategy using various thresholds. These columns are further subdivided into two, 
the first containing X^n=i dimV™ (i.e., the total number of degrees of freedom from 
all meshes over time) which we denote DOF and the second containing CPU time 
(sees) for all model problems (|5.1|)-(5.3|. 





Uniform 


Adaptive 






£ = 0.65 


5 = 0.70 


( = 0.75 


tol 


DOF's 


CPU 


DOF's 


CPU 


DOF's 


CPU 


DOF's 


CPU 


0.573 


232,290 


3 


24,080 


4 


22,792 


5 


22,240 


4 


0.295 


3,489,090 


49 


42,042 


8 


39,414 


8 


38,630 


6 


0.149 


54,097,020 


598 


82,172 


15 


77,932 


15 


76,452 


16 


0.0625 


OOM 


OOM 


206,709 


39 


195,810 


37 


191,650 


37 



Table 1 . Explicit timestep control with various spatial maximum 



strategy thresholds for Problem (5.1 1. The adaptive method clearly 



saves DOF and CPU time over the uniform method. 





Uniform 


Adaptive 






f = 0.65 


5 = 0.7 


£ = 0.75 


tol 


DOF's 


CPU 


DOF's 


CPU 


DOF's 


CPU 


DOF's 


CPU 


0.296 


3,489,090 


47 


12,092 


5 


11,430 


5 


11,498 


5 


0.21 


13,940,289 


196 


17,038 


7 


10,140 


8 


16,201 


7 


0.104 


54,097,020 


602 


106,188 


32 


100,058 


29 


22,597 


10 


0.03125 


OOM 


OOM 


513,694 


120 


460,637 


118 


449,568 


115 



Table 2. Explicit timestep control with various spatial maxi- 



mum strategy thresholds for spatial-error dominant Problem (5.3) 
Adaptivity saves DOF and CPU. 
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Uniform 


Adaptive 






£ = 0.7 


£ = 0.75 


tol 


DOF's 


CPU 


DOF's 


CPU 


DOF's 


CPU 


1.000 


925,809 


12 


159,070 


43 


127,610 


58 


0.569 


3,489,090 


49 


237,960 


142 


204,376 


180 


0.295 


54,097,020 


605 


471,733 


755 


471,542 


920 


0.149 


OOM 


OOM 


940,618 


1410 


940,138 


1850 



Table 3. Implicit timestep control with various spatial maxi- 



mum strategy thresholds for spatial-error dominant Problem (5.2). 
Adaptivity saves DOF (even better than explicit control) but the 
CPU time grows very quickly due to overhead. 





Uniform 


Adaptive 






£ = 0.7 


£ = 0.75 


tol 


DOF's 


CPU 


DOF's 


CPU 


DOF's 


CPU 


1.000 


925,809 


12 


135,788 


5 


127,004 


4 


0.569 


3,489,090 


49 


198,628 


7 


194,311 


8 


0.295 


54,097,026 


605 


397,716 


15 


395,876 


16 


0.149 


OOM 


OOM 


2,177,666 


79 


2,079,081 


76 



Table 4. Explicit timestep control with various spatial maximum 



strategy thresholds for time-error dominant Problem (5.2 1 



6.7. Remark (implicit timestep control on fast oscillating solutions). We take note 



of the CPU times from the results for Problem (5.2 1. These show that implicit 
timestep control is undesirable for fast oscillating functions. This is because the 
timestep searching becomes computationally inefficient. Numerical simulations for 
an explicit timestep control strategy is given in Table |4j This algorithm is described 
in detail in the ALBERTA manual SS05, §1.5.4] The results show although for a 
method with low tolerance we use more degrees of freedom we make a substantial 
gain on the CPU time. 

We then fix a value of £ and compare an adaptive strategy with uniform for a 
single value of tol. This is to illustrate how the number of degrees of freedom of the 
mesh change over time, and how the implicit timestep control affects the timestep 
size for all test problems in Figures [6| 

6.8. Incompatible data singular solution. We close the paper by testing the 
adaptive algorithm on an example with incompatible initial and boundary condi- 
tions, which is the type of situation where adaptivity is really needed in practise. 



Consider problem (2.6| with Q = (0, 1) x (0, 1), / = and uo = 1. The initial con- 
ditions are thus incompatible with the homogeneous Dirichlet boundary valid for 
all positive times. The exact solution u, though singular at all points of df2 x {0}, 
can be readily evaluated "by hand" and may be represented in terms of Fourier 
series of the Laplacian's eigenvalues. Namely, we have 



oo 

u(x, t) = C„ lyn exp(— (m 2 + n 2 ) ir 2 t) sin(m7rai) sin(n7TX2), for t > 0, (6.2) 

m,n— 1 



where the constant C m>n is given by 
4 

C m ,n = J (1 ~~ C0S ( TO7r ) — C0S ( n7r ) + C0S ( n7r ) C0S (WWI")) . (6-3) 



Since the solution (6.2) is an infinite Fourier series it cannot be computed exactly, 



but its rapid decay allows to truncate early with machine-epsilon precision. 
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In order to generate a reference tolerance, which is common for the uniform and 
the adaptive scheme we couple h = 0.05r and run the uniform refinement code. We 
use then the error computed as a tolerance for the adaptive scheme, results of this 
are shown in Figure [7j In Figure [8] we visualise the adapted FE mesh for Problem 



(6.2) at various times. 
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Appendix A. Coarsening error preindicator implementation 

We describe next a practical implementation of the coarsening error preindicator 
(we use this term to emphasise the fact that this indicator can be computed a priori, 
as opposed to the other indicators involved in the adaptive strategy) . Since we used 
ALBERTA for our computations, this section relies substantially on the principles 
described in the manual |SS05j . We briefly describe these principles in the next 
paragraph, in order to expose the main idea behind the coarsening preindicator. 

A.l. Refinement, coarsening and interpolation in ALBERTA. Mathemati- 
cally, a simplicial mesh (or partition, or triangulation) is a set of disjoint open 
simplexes, the union of the closure of which is Q. A mesh into a new mesh is 
refined by bisecting a subset of its simplexes, following a special procedure which 
ensures mesh conformity (e.g., no hanging nodes) and does not deteriorate shape- 
regularity (on fully fitted polygonal domains). A mesh is thus represented as a 
binary tree, where each node represents a simplex. The children of each simplex 
are thus the 2 subsimplexes obtained by bisection. Hence, from a coding view-point, 
refinement means growing the binary tree. 

The inverse of refinement is coarsening. Thus coarsening a mesh in ALBERTA 
consists in removing pairs of sibling simplexes (both marked for coarsening) and 
produces the new — coarsened — mesh where the pairs of siblings are replaced by 
their parent. 

The coarsening preindicator is a real number defined on each simplex, of the tri- 
angulation to be coarsened. This estimator can in fact be precomputed with respect 
to coarsening. This is in contrast with usual a posteriori error estimators which 
can be postcomputed only (i.e., after the discrete solution has been computed). To 
clarify this point, let us focus on the particular situation of interest. Let U n ~~ x be 
the solution from the previous timestep; U n ^ 1 € V™ -1 , the finite element space 
with respect to mesh , ( 7 n ^ 1 . The error due to coarsening appears in the term 

This term is nonzero only when simplexes are coarsened. 

Furthermore, we assume that the new mesh ST 71 is a refinement of =^q™, which is 
a coarsening of the old mesh 3" a ^ 1 : 

Ofn-l coarsen^ ^ n refine^ _ _ _ refine^ ^ n ^ ^\ 

If Aq is the Lagrange interpolator onto the finite element space Vq, relative to 
the new coarse mesh <5^ n , it is not very difficult to predict /IqC/™ -1 without actually 
computing it. Therefore this term can be predicted from (a) the simplexes of 3~ n ~ x 
marked for coarsening which leads to 3?^" and (b) the values of [7 n_1 . 

Note that since S?^ 1 is subsequently refined but not coarsened to produce ST n , as 



depicted in (A.2|, then the additional coarsening error will be zero. Namely, if A 1 
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denotes the Lagrange interpolant onto V™, the finite element space over 3? n , which 
is a refinement of ^ ", then A n U n - x = A^U n - 1 , and thus 

Tjn-l _ A n un -1 = jjn-l _ J^jju-l ^gj 

The coarsening strategy therefore consists in choosing a subset of simplexes of 
^ n_1 which minimises term — t/" - 1 1 1 before producing the new coarse 

mesh % n . 

The rest of this section describes how C/™ -1 — t1q£/™ -1 can be precomputed. 

A. 2. Notation. Let K be an element of the new coarse mesh 3? n resulting from 
the coarsening of its two children which we denote by K ± . (Note that K + and 
K~ correspond to child [0] and child [1] of K in the ALBERTA manual jSS05j .) 
Define the fine space 

Y := {<P\ K : $ e V"" 1 } . (A.4) 

Likewise define the coarse space X to be the local finite element space, i.e., 

X :={<%- : <Z>eV£}; (A.5) 

simply put we just have X = P p . We introduce also the fine spaces Y 1 * 1 , defined like 
Y, but restricting functions over K ± , respectively (so functions in Y* are in fact 
the same as X = P p , albeit with different domains). 

Denote by {x , . . . , x^} and {x J , . . . , x ^ } the set of Lagrange degrees of freedom 
on the simplex K and its children K , respectively. We indicate with {ir°, . . . , 7r L } 
and {tt±, . . . ,tt±} the corresponding Lagrange polynomial bases of X and Y*, re- 
spectively, whereby 

n i (x j )=7r i ± (xf) = 6 i j . (A.6) 

For short we will write these bases as column vectors 7r = (tt , . . . ,7r L ) T , etc. We 
also define the (local) coarse-on-fine matrixes by 

A± := ( W (x±) . . . = (A*f)) id=0t .., tL ■ (A.7) 

These matrixes are closely related to ALBERTA's refine-interpolation matrix [SS05 , 
matrix A (1.5) in §1.4.4 ]. 

A. 3. Proposition (coarse-on-fine matrix properties). The matrixes A + and A arc 
independent of K, K + , K~ and 

7r|^ ± =A ± 7 r ± . (A.8) 

Proof Fix i = 0, . . . ,L. Because it 1 is a polynomial and {ir^ , . . . ,7! + } is a polyno- 
mial basis, it follows that 

L 

^ = 5>K> (A.9) 

j=0 



for some vector (a , . . . ,a l L ). Applying n l to x^ , and recalling (A.6), we obtain 



rr l (x+), (A. 10) 



and hence 



[A+TT+f. (A.ll) 

□ 



A.4. Example (quadratic elements in 2 dimensions). To make the discussion more 
accessible, we will illustrate it as we go with the concrete situation where p = 2 
(quadratic elements) and d = 2. Following the ALBERTA conventions the relation 
between the coarse and fine triangles is given by the following diagram. 
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coarsen(X + , K 



® ® 



K+ K K~ refined K+ ' K~ 

© © ® ® © ® ® < 

In this case, the coarse-on-fine matrixes are computed as follows: 



3/8 -1/8 
-1/8 -1/8 

1/2 

1/2 1 

3/4 1/4 



-1/8 -1/8 
-1/8 3/8 

1/2 
1/2 

1/4 3/4 



(A.12) 



A. 5. Degrees of freedom and global— local relations. Denote by U the generic 
finite element function in the old space V™ -1 and let V := AqU. Then we have 

U = u T «P and V = v T <P, (A. 13) 

where \P = (if -0 , . . . ,!f rAr ) T and # = (<P°, . . . ,<£ M ) T , are the columns of nodal La- 
grange piecewise polynomial bases of V 71-1 and Vq , respectively, and u and v are 
the corresponding vectors of DOF values. 

There are L + l degrees of freedom (DOF) per simplex, e.g., L — 5 for p = 2 = d. 
The simplex K in 2?^ comes with a local-to-global index relation g — g K ° : [0 : L] — > 
[0 : M) whereby 



K 



7T* Vj=0, 



,L. 



(A.14) 



It follows that the finite element function V is locally represented on K by 

L 



i=0 



Similarly we have g = g^± '■ [0 : L] [0 : N] such that 



Y ± := U\ K± = ^w s ± 0) 7ri =: y ±T 7r±. 

3=0 

The relation between the DOF coefficients u and v will be described next. 



(A.16) 



A. 6. Local fine— coarse DOF relations. Some degrees of freedom — that is those 
depicted in yellow or bright — are removed during coarsening. The others, which 
are kept, have their local index change. This information is fully encoded in the 
fine-to-coarse index maps : D — > C where 

D ± :={j = 0,...,L: xf € {x , . . . ,x L }} . (A.17) 

and 

^ :=c ± (D ± )C[0:L}. (A.18) 

A basic property of the fine-to-coarse maps is that 

C+UC- = [0: L], (A.19) 

but C + and C~ need not be disjoint (in fact, for conforming methods these are never 
disjoint). The fine-to-coarse maps c ± are injective and we denote their inverses, the 
coarse-to-fine maps, by d : C — > D . 
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In the example above, p = 2 = d, the fine-to-coarse maps c ± : D ± — > [0:5], 
satisfy D + = D~ = {0, 1,2,5} (though D + and D~ do not generally coincide, as 
seen for p = 3, d = 2, e.g.) and evaluated by the schedule 

.7 = 1 2 3 4 5, 
c+(j) = 2 5 - - 4, (A.20) 
c~(j)= 1 2 5 - - 3. 

It follows that C+ = {0, 2, 4, 5} and C~ = {1, 2, 3, 5} and 

i= 1 2 3 4 5, 
=1-0-52, (A.21) 
=-015-2. 

A. 7. Remark (redundancy of the coarse-to-fine maps). The coarse-to-fine maps 
and their inverses d^ are partially redundant with A ± . Namely, if j 6 Z^, then 

j = d ± (i) and i = c (j), for some i = 0, . . . , L. By definition of c it follows that 
= £Ej. Therefore 

[A^tt*^) (A.22) 

We have thus proved the following result that will be used to compress A* in the 
sequel. 

A. 8. Proposition (redundant coarse-on-fine columns). If j 6 D ± , then A^'sj-th 
column is described by 

[A% fc =<5 c fc ±w . (A.23) 

A. 9. Precomputing the coarsening error. The coarsening error is the differ- 
ence between U, to which we have access via u, and its interpolation on the lo- 
cally coarser mesh V, to which we have no direct access. Working locally at the 
coarsening-marked element K + (and similarly for K^) 1 all we need is to compute 
^1^+ and subtract it from U\ K+ . 

Recalling that in ALBERTA V = A^U is built by simply "dropping" the coeffi- 
cients of the DOF removed by coarsening we have 

u g+ (d+(i))K l + ^2 U g(d - (l)) 1T\ (A.24) 

ie 

= tf*ec+ (A25) 

*a-(«J-W) =y d -(i) otherwise . 

(Note that the vector y is the same for the two siblings K ± and needs to be 
calculated only once.) Following the example with p = 2 = d, we see that 

y = {vtiVa >yt>yh ,vtiviY 
= {yf ,Vo ,vi ,y$ , y£,vi 

To conclude we rewrite the coarse basis, 7r, in terms of the fine one, 7r + , using 
Proposition |A.3| as follows: 

V\k+ = y \k + = y T = y T A + 7r+. (A.27) 

Thus the coarsening error on K + is calculated as 

L 

P - V]\ K+ = y+V - yTA+7r + = 7r + T(y+ - A +T y) = £(y+ - yT[A+] 3 ) tt+ 

3=0 

(A.28) 




6C+ i£C-\C+ 



(A.26) 
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Recalling Proposition A. 8 if j <E D + we have 

L 

Y1A+] . = Vk&w = Vo+U) = Vf, (A.29) 

j=0 

and thus the coefficient for is 0, and it needs not be calculated. Proceeding 
similarly on K~ we may summarise the findings as follows. 

A. 10. Theorem (coarsening error calculation). Let U € V™ -1 with the notation 
of §A.5[ to calculate the coarsening error that would result from coarsening the 
elements K + ,K~ G into K G ,9~ n 



(1) calculate y following \A.25) using the coarse-to-fine map d + defined in [A.6 

(2) obtain the error using 

[u-a%u\\ k+ = Yl {vf-yWi)*- 

jE[0:LWD+ 

(A.30) 

[u-a%u]\ k -= (MT-y^'W**- 

je[Q:L]^D- 

A.ll. Remark. Note that the j-th coefficient of the coarsening error's local DOF 
vector is zero when j G D , respectively. So the calculation needs to be carried 
out only for those j £ D ± . 

Also, the coefficients for the DOF that are common to K + and K~ must be 
equal, so they can be in fact computed once. 

For example in the case of quadratic elements in d — 2 we have 
Y + -Y\ K+ =nl(y+ 



Y~-Y\k- =^-1% 



I 2/4 



3 + 


+ 8% 








1 + 

+ 8 yt 


+ 8% 




2% 


1 


i + 
+ - 8 yt 




1 


2% 


V 


i + 

+ - & vt 


3 









(A.31) 



A. 12. Coarsening error algorithm. As seen in |A.9| the information needed 
for the coarsening error computation for Lagrange finite elements of d egree p in 
dimension d, is contained in the coarse-on-fine matrixes A ± defined by (A.7) and 
the fine-to-coarse maps, d^, and their domains C ± defined in A.6 This information 
is independent of the particular pair of simplex siblings K ± and their parent K and 
can be included in the code via given index permutations and efficient matrix- vector 
multiplication. 

With this information at hand and the notation previously introduced in this 
section, we formulate an ALBERTA-implementable algorithm to precompute the 
coarsening error on all seimplexes. 

Coarsening Preindicator. 

Require: ([/ = u T *,V, ST) 
Ensure: 7 = (*y K : K G 
for all K G & do 

if childorder(lf) = C0 then 



2 The element information in ALBERTA is quite local and to determine whether an element 
is left or right child is not trivial. In ALBERTA 1.2 this can be done utilising EL->index which 
provides a global indexing of elements. Testing the EL_INFO->parent->child [0] ->index against 
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:= c+ d := c , A := 



D := D+, D' := D~ 
else 

D := D~, D' := D+, c := c~, d := c 
end if 

K' := sibling if 

initialise two local DOF vectors y and r 
for all j e D do 

Vc(j) = U 9kU) 

end for 

for all j e D' do 

ycu) = u 9k>u) 

end for 

for all j <£ D U D' do 

r j = u 9K(j) -y T I A ]j 

end for 

7k = 

for alii £ D U D' do 
for all j^DUfl'do 

IK = lK + r l r J (<P l: <Pj) K 
end for 
end for 
end for 



A. 13. Coarsening preindicator matrixes. To close, we provide here the in- 
formation needed to implement Algorithm |A.12| for Lagrange piecewise P p finite 
elements in dimension d — 2. (For dimension 3 the situation is complicated by the 
"types" of tetrahedrons, whereby the matrixes A ± and the maps c may depend 
on the type and is not covered in this appendix.) 



A. 14. P 1 elements. The coarse-on-fine matrixes (omitting entries for clarity) 
are given by 





1 1/2" 




r i/2 


A+ = 


1/2 


A" = 


1 1/2 




1 




1 



(A.32) 

the fine-to-coarse maps and the coarse-to-fine maps are respectively given by 

i = 1 2, i = 1 2, 

c+(i) = 2 -, and d+(i) =1-0, (A.33) 
c (i) =12-, d (i) =-01. 



A. 15. P 2 elements. See the worked example in §A[ 



EL->index gives the correct child order of K. In ALBERTA 2.0 EL->index is unavailable so we 
check the global index of DOF for both parent and children. 
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A. 16. P 3 elements. The coarse-on-fine matrixes are given by 



A+ = 



and 



-1/16 5/16 
-1/16 1/16 



9/16 15/16 1 
9/16 -5/16 



-1/16 
-1/16 



9/16 
9/16 



the fine-to-coarse maps 

c+(i) = 
c-(i) = 

and the coarse-to-fine maps 



1/16 
1/16 

-1/4 

1/2 

1/2 
-1/4 
-1/16 
-1/16 

1/2 

1 2 

- 

2 - 



(i) 



1 

1 - 
- 



1/16 
1/16 

-1/4 
1/2 
1/2 
-1/4 
-1/16 
-1/16 
1/2 



1/16 
5/16 



-5/16 
15/16 



-1/16" 
1/16 

-1/8 



3/8 
3/16 
-3/16 
3/4 



-1/16" 

1/16 

-1/8 
-3/16 

3/16 

3/8 

3/4 



7 8 9, 

5 6-, 

3 4-. 

7 8 9, 

4 5, 

- 5 4. 



(A.36) 



(A.37) 
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A. 17. P 4 elements. The coarse-on-fine matrixes are given by 



00 00 
CM CM 



in in 
I 



co in 
I 



CO 00 
I 



tO CM 3f CM 2 «3 
i — I CO cP CO i — I 



« S co °° °° 

00 ^ oo co~ m ="? 



00 GO 
CM CM 



in co 
I 



00 CO 
CM CM 



in in 
I 



. '20' 

co 



CO CM 3 CM CO 2 

i — i co -cr^ co i — i 

I 



co oo 

CM CM 



in in 
I 



CO I 



CO 

I 



H CO CO 

I 



00 CM 
CM 2 
— [ ' 1 



CM 
CO CO 

• in in 
co co 



to cs n cs oi 
I 



~j -tf CM 
CO C° « 



in in 
^ ' — i 



CO 00 
CM CM 



CM 

co 



00 CO 
CM CM 



in in 



CD CM CM CO „ 

H CO tO M H © 



H M Ol M ffi 

I 



CO M 
CM N 



in "5 
co | 



2 00 CM CMOO 2 « S CM CO CO ^ 



in in 
co co 
I 



in in 
I 



■CO H H CO 



i — I CO CO 

I 



+ 
< 



The fine-to-coarse maps are given by 

i= 1 2 3 4 5 6 

c+(i) =20 10- 9 - - 

c~(i) = 1 2 10 - 14 - - 

and the coarse-to-fine maps by 

i= 1 2 3 4 5 6 

d+(i) = 1 - - - - 9 

d (i)= - 1 9 10 11 - 



7 8 9 10 11 12 13 14, 

14 - 6 7 8 - - 12, 
11 - 3 4 5 - - 13. 



7 8 9 10 11 12 13 14, 
10 11 4 2 - 14 - 7, 
- - - 2 7 - 14 4. 
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FIGURE 1. Numerical Results for P roble m (5.1) with P 1 and h 
h(i) = 2~ 1 / 2 , i = 4, . . . , 9 (details in jj5j}. 
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EOC[E 



EOC[e ] 



5 
4 
3 
2 
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(a) Mesh-size is h and timestep r = 0.1 h. On top we plot the EOC's of the single 
cumulative indicators E and 0. Below we plot their logs. Both indicators have EOC — > 1, 
but the cumulative time error indicator 6> m is dominant. The estimator is reliable and 
sharp, but not asymptotically exact and results in EI 3> 1. 




■ 0' ■ 0' ■ 0' ■ 

0.5 1 0.5 1 0.5 1 0.5 1 




0.5 1 0.5 1 0.5 1 0.5 1 



(b) Timestep is t = O.lh 2 . This choice leads to EOC[<9 m ] -> 2 and EOC[E m ] 1, i.e., 
the time indicator m is of higher order than the spatial indicator E m which leads the 
estimator's order. Thus we obtain asymptotic exactness EI — > 1, as expected from ZZ 
estimators for p = 1. 
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(a) Timestep r = 0.1 h 2 . The cumulative time error indicator <9 m is dominant with 
EOC[6> m ] -> 2, but EI > 1. 




(b) Timestep is r = 0.1 h 3 , with In the bottom set of results the spatial is dominant 
(EOC ~ 2) showing the estimator is sharp and reliable for higher order polynomials as 
well, and close to asymptotically exact (EI just smaller than 1). 
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Figure 3. Numerical Results for (5.1 ) with P 3 elements for mesh- 
sizes h(i) — 2~ 1 / 2 , i = 2, .. .,6.. We compute the same quantities 
as in Figures [T] and [2j 



E0C X ,m w $ Eoc[E m ; 

5r ■ ■ : 5 



Eoqe ] 




0.5 1 0.5 1 0.5 1 0.5 1 




0.5 1 0.5 1 0.5 1 0.5 1 



(a) Timestep is r = 0.1 ft 3 . Again, the time indicator is dominant and EOC[<9 m ] — > 3, 
but EI > 1. 



EOC[E 



Eoqe ] 



0.5 1 0.5 1 0.5 1 





0.5 1 0.5 1 0.5 1 0.5 1 



(b) Timestep is r = 0.1 ft 4 . The elliptic error is dominant (EOC[i? m ] — > 3) and the 
estimator is sharp and reliable with very good EI. 
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FIGURE 4. Results for (|5.1[) with P 4 elements and h(i) = 2~ 1 / 2 , 



i = 2, . . . , 6. We compute the same time accumulation quantities 
as in Figures [T] |3] 



E0C [" ell L ! « l , -; H;,] EOC[E m 

5 5 



2L 
1 



E0C[9 ] 





0.5 1 0.5 1 0.5 1 0.5 1 




0.5 1 0.5 1 0. 



(a) Mesh -size is t — 0.1 . Again, the time indicator is dominant with order EOC[@tj 
4) and a quite good EI in this case. 



EOC[llell 




(b) Mesh-size is r = O.l/i . The spatial error is dominant and EOC[£? m ] — > 4. EfTectivity 
index improves slightly over previous case. 
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Figure 5. For each m = 1,...,N we plot values and EOC's 
of two alternative time indicators (Y]™—i Tn&^J (above) and 

(Sr=i Tn^n) 1 ^ 2 (below) and the alternative mesh-change indicator 
E™=i T nln (above-right). All quantities are plotted against time. 
We took a uniform timestep r = 0.1 h and mesh-size h = 2~ % , 
i = 4, ... ,9. The numerical results show (1) that the two time 
indicators are equivalent in order, as expected, and (2) that the 
term X^=i T »7n * s indeed a higher order term and can be safely 
ignored in most practical schemes. The indicators 9 n have a better 
effectivity index. 



2.5 r 2.5 




0.5 1 0.5 1 0.5 1 0.5 1 
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Figure 6. Adaptive (green) against uniform (red) degrees of free- 
dom and timestep sizes. In each pair of graphs we plot the (log of) 
the DOF against time on the left, and the timestep against time 
on the right. 





(a) Implicit timestep control for Problem 
l |5.1| . The explicit timestep control yields 
the same results (but is much more CPU 
efficient), thus it is not shown. 



(b) Implicit timestep control for Prob- 
lem l ]5.2| . where the spatial error domi- 
nates. The explicit timestep control yields 
the same meshes and time-steps, thus not 
shown. 



UMUUUUUUULUUUUUUl 



UIMMJUUIMMMJUMI 



(c) Explicit timestep control for Problem 
l |5.3| , where the time discretisation error 
dominates. Interesting when compared 
with Figure [7(d)] 



(d) Implicit timestep control for Problem 
fl5.3| . Comparing with Figure [7 (c) | shows 
that the implicit timestep control yields 
more efficient timestep and meshes, but 
at a much higher CPU cost (cf. Tables [3] 
and [If. 
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Figure 7. Implicit timestep control for Problem (6.2 1. 




Figure 8. The adaptive scheme for (6.2 1 using implicit timestep control 






mm 



(a) Mesh at time (b) Mesh at time (c) Mesh at time (d) Mesh at time 
t n = 0.007544 with t„ = 0.033302 with t n = 0.127492 with t n = 0.393893 with 
dim(V") = 894, 677 dim(V") = 98, 773 dim(V n ) = 18, 613 dim(V") = 3, 525 



